tidymodels를 이용한 시계열 모델링
modeltime 패키지는 tidymodels와 연동이 가능한 시계열 모델링 관련 패키지이다. tidymodels에도 auto arima, ma 모형 같은 간단한 시계열 모델이 있지만 다른 머신러닝 모델에 비해서 쓸 수 있는 모델이 제한적이다. 이러한 단점을 modeltime 패키지가 해결해준다.
modeltime 시계열 모델링에 특화된 패키지로 시계열에 특화된 모델링과 전처리 관련 함수가 내장되어 있고, tidymodels의 워크플로우를 거의 그대로 이용할 수 있어서 향후 tidymodels가 R의 대표적인 머신러닝 패키지가 된다면 시계열 파트에서는 modeltime 패키지가 주축이 될 것 같다.
modeltime 패키지는 하나의 단일 패키지가 아니라 다양한 머신러닝 패키지와 연동해서 하나의 시계열 생태계를 구축하고 있다. 대표적인 패키지는 다음과 같다.
Modeltime : 시계열 머신러닝 관련 패키지
Modeltime H2O : H2O의 autoML과 연동이 가능함
Modeltime GluonTS : 시계열 관련 딥러닝 패키지
Modeltime Ensemble : Modeltime 관련 앙상블 패키지
Modeltime Resample : Backtesting 관련 패키지
Timetk : feature engineering, Data wrangling, time series visualization
Modeltime 패키지는 시계열 모델링에 맞는 고유한 워크플로우가 존재하며, tidymodels와 마찬가지로 이 워크플로우 그대로 진행을 해야만 에러 없이 동작한다. 사실 실제 모델링할 때와 동일한 매커니즘이어서 어색함 없이 따라할 수 있다.
workflow의 순서는 다음과 같다.
create modeltime table
다양한 모델을 training data에 fitting하고 modeltime table에 저장하는 단계
모델이 재대로 fitting 되었는지 확인 및 이후 워크플로우에 맞는 table 구성
calibrate
modeltime table에 저장된 모델을 test 데이터에 fitting하는 단계
test data에 fitting함으로써 모델의 성능 파악 가능
Refit
training data, test data를 합친 original data에 refit하는 단계
refit한 모델을 이용해서 forecast를 진행
Modeltime 패키지는 다양한 시계열 모델을 지원하는데 대표적인 모델은 다음과 같다.
Modeltime 패키지의 한 가지 단점은 패키지를 만든 저자가 business-science라는 course를 운영 중인 것 같은데 modeltime 관련 세부 자료는 유료 course에서 공개를 하는 것 같다… 하지만 R CRAN에도 등록이 되어있고 튜토리얼도 제공되기 때문에 앞으로 사용하는 사람이 많아지면 관련 자료도 많아질 것 같다.
데이터는 데이콘에서 제공하는 발전량 데이터 일부를 뽑아서 활용했다.
library(tidyverse)
library(tidymodels)
library(lubridate)
library(data.table)
library(skimr)
library(timetk)
library(modeltime)
library(gt)
library(timetk)
library(tidyquant)
library(sknifedatar)
# Visualization
library(ggthemes)
library(ggsci)
library(viridis)
library(ggExtra)
theme_set(theme_bw())
file_path <- getwd()
files <- list.files(file_path)
files
[1] "modeltime-with-tidymodels.html"
[2] "modeltime-with-tidymodels.Rmd"
[3] "modeltime-with-tidymodels_files"
[4] "rdata.csv"
head(rdata)
time 기온 강수 풍속 풍향 습도 증기압 이슬점 기압
1: 2015-01-01 01:00:00 -4.4 0.00 5.4 340 47 2.1 -14.0 1020.3
2: 2015-01-01 02:00:00 -4.6 0.00 4.9 340 50 2.2 -13.4 1020.3
3: 2015-01-01 03:00:00 -4.7 0.05 6.2 320 50 2.2 -13.5 1020.7
4: 2015-01-01 04:00:00 -5.0 0.00 5.0 320 56 2.4 -12.4 1020.6
5: 2015-01-01 05:00:00 -5.0 0.00 5.5 320 52 2.2 -13.3 1020.4
6: 2015-01-01 06:00:00 -5.3 0.00 4.0 340 58 2.4 -12.2 1020.7
해면기압 일조 일사 전운량 시정 지면온도 floating warehouse dangjin
1: 1024.0 NA 0 6 1500 -4.4 NA 0 0
2: 1024.0 NA 0 NA 1750 -4.6 NA 0 0
3: 1024.5 NA 0 6 2000 -4.7 NA 0 0
4: 1024.4 NA 0 6 2000 -5.0 NA 0 0
5: 1024.2 NA 0 6 2000 -5.0 NA 0 0
6: 1024.5 NA 0 6 2000 -5.3 NA 0 0
ulsan hour
1: 0 1
2: 0 2
3: 0 3
4: 0 4
5: 0 5
6: 0 6
glimpse(rdata)
Rows: 55,392
Columns: 20
$ time <dttm> 2015-01-01 01:00:00, 2015-01-01 02:00:00, 2015-01~
$ 기온 <dbl> -4.4, -4.6, -4.7, -5.0, -5.0, -5.3, -5.7, -5.7, -5.4~
$ 강수 <dbl> 0.00, 0.00, 0.05, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00~
$ 풍속 <dbl> 5.4, 4.9, 6.2, 5.0, 5.5, 4.0, 5.0, 4.6, 6.7, 6.5, 6.~
$ 풍향 <dbl> 340, 340, 320, 320, 320, 340, 340, 340, 320, 320, 32~
$ 습도 <dbl> 47, 50, 50, 56, 52, 58, 58, 56, 57, 59, 60, 62, 58, ~
$ 증기압 <dbl> 2.1, 2.2, 2.2, 2.4, 2.2, 2.4, 2.3, 2.3, 2.3, 2.5, 2.6~
$ 이슬점 <dbl> -14.0, -13.4, -13.5, -12.4, -13.3, -12.2, -12.6, -13.~
$ 기압 <dbl> 1020.3, 1020.3, 1020.7, 1020.6, 1020.4, 1020.7, 1021~
$ 해면기압 <dbl> 1024.0, 1024.0, 1024.5, 1024.4, 1024.2, 1024.5, 1025.1~
$ 일조 <dbl> NA, NA, NA, NA, NA, NA, NA, 0.0, 0.0, 0.6, 0.9, 0.5,~
$ 일사 <dbl> 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.000, 0.0~
$ 전운량 <int> 6, NA, 6, 6, 6, 6, 6, 3, 6, 6, 7, 8, 8, 8, 8, 6, 6, 8~
$ 시정 <dbl> 1500.000, 1750.000, 2000.000, 2000.000, 2000.000, 20~
$ 지면온도 <dbl> -4.4, -4.6, -4.7, -5.0, -5.0, -5.3, -5.7, -5.7, -5.4, ~
$ floating <dbl> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA~
$ warehouse <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 12, 150, 195, 270, 254, 18~
$ dangjin <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 8, 227, 263, 356, 333, 225~
$ ulsan <int> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,~
$ hour <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15,~
skim(rdata)
| Name | rdata |
| Number of rows | 55392 |
| Number of columns | 20 |
| Key | NULL |
| _______________________ | |
| Column type frequency: | |
| numeric | 19 |
| POSIXct | 1 |
| ________________________ | |
| Group variables | None |
Variable type: numeric
| skim_variable | n_missing | complete_rate | mean | sd | p0 | p25 | p50 | p75 | p100 | hist |
|---|---|---|---|---|---|---|---|---|---|---|
| 기온 | 0 | 1.00 | 12.21 | 10.36 | -19.3 | 3.60 | 12.60 | 21.10 | 36.80 | ▁▅▇▇▂ |
| 강수 | 0 | 1.00 | 0.12 | 1.05 | 0.0 | 0.00 | 0.00 | 0.00 | 102.70 | ▇▁▁▁▁ |
| 풍속 | 0 | 1.00 | 1.96 | 1.56 | 0.0 | 0.70 | 1.60 | 2.90 | 11.70 | ▇▃▁▁▁ |
| 풍향 | 0 | 1.00 | 164.29 | 129.50 | 0.0 | 20.00 | 180.00 | 290.00 | 360.00 | ▇▁▃▂▆ |
| 습도 | 0 | 1.00 | 75.34 | 20.14 | 10.0 | 61.00 | 79.00 | 94.00 | 100.00 | ▁▂▃▅▇ |
| 증기압 | 82 | 1.00 | 12.93 | 8.92 | 1.0 | 5.40 | 10.20 | 19.20 | 42.90 | ▇▅▃▂▁ |
| 이슬점 | 85 | 1.00 | 7.37 | 11.04 | -22.4 | -1.70 | 7.30 | 16.90 | 30.20 | ▁▆▇▇▅ |
| 기압 | 79 | 1.00 | 1014.05 | 8.36 | 983.6 | 1007.40 | 1014.30 | 1020.60 | 1036.30 | ▁▃▇▇▂ |
| 해면기압 | 78 | 1.00 | 1017.37 | 8.48 | 986.6 | 1010.60 | 1017.60 | 1024.00 | 1039.70 | ▁▃▇▇▂ |
| 일조 | 25453 | 0.54 | 0.52 | 0.45 | 0.0 | 0.00 | 0.60 | 1.00 | 1.00 | ▇▁▁▁▇ |
| 일사 | 0 | 1.00 | 0.54 | 0.83 | 0.0 | 0.00 | 0.01 | 0.89 | 4.85 | ▇▂▁▁▁ |
| 전운량 | 12427 | 0.78 | 5.23 | 3.84 | 0.0 | 1.00 | 6.00 | 9.00 | 10.00 | ▇▂▃▅▇ |
| 시정 | 0 | 1.00 | 1754.20 | 1024.17 | 3.0 | 1100.00 | 1800.00 | 2029.00 | 6454.00 | ▅▇▁▁▁ |
| 지면온도 | 0 | 1.00 | 12.21 | 10.36 | -19.3 | 3.60 | 12.60 | 21.10 | 36.80 | ▁▅▇▇▂ |
| floating | 28344 | 0.49 | 122.42 | 192.34 | 0.0 | 0.00 | 0.00 | 192.25 | 753.00 | ▇▁▁▁▁ |
| warehouse | 2040 | 0.96 | 95.81 | 150.72 | 0.0 | 0.00 | 0.00 | 152.00 | 595.00 | ▇▁▁▁▁ |
| dangjin | 2040 | 0.96 | 140.11 | 221.67 | 0.0 | 0.00 | 0.00 | 227.00 | 881.00 | ▇▁▁▁▁ |
| ulsan | 3528 | 0.94 | 66.28 | 104.17 | 0.0 | 0.00 | 0.00 | 104.00 | 448.00 | ▇▁▁▁▁ |
| hour | 0 | 1.00 | 11.50 | 6.92 | 0.0 | 5.75 | 11.50 | 17.25 | 23.00 | ▇▇▆▇▇ |
Variable type: POSIXct
| skim_variable | n_missing | complete_rate | min | max | median | n_unique |
|---|---|---|---|---|---|---|
| time | 0 | 1 | 2015-01-01 01:00:00 | 2021-04-27 | 2018-02-28 00:30:00 | 55392 |
rdata %>%
select(-hour) %>%
mutate(time = ymd_hms(time)) %>%
filter(between(time, ymd('2018-03-01'), ymd('2021-01-31'))) -> rdata
rdata %>% glimpse()
Rows: 25,609
Columns: 19
$ time <dttm> 2018-03-01 00:00:00, 2018-03-01 01:00:00, 2018-03~
$ 기온 <dbl> 3.1, 2.8, 2.6, 2.0, 2.2, 4.1, 3.5, 2.2, 1.0, 0.3, 0.~
$ 강수 <dbl> 0.50, 0.00, 0.00, 0.00, 0.00, 0.00, 1.80, 0.00, 0.00~
$ 풍속 <dbl> 3.6, 0.7, 3.2, 1.9, 2.1, 4.4, 7.9, 6.4, 7.7, 8.9, 7.~
$ 풍향 <dbl> 340, 140, 320, 230, 180, 270, 320, 290, 320, 320, 32~
$ 습도 <dbl> 96, 97, 95, 97, 97, 97, 93, 86, 82, 71, 63, 58, 60, ~
$ 증기압 <dbl> 7.3, 7.2, 7.0, 6.8, 6.9, 7.9, 7.3, 6.1, 5.4, 4.4, 4.0~
$ 이슬점 <dbl> 2.5, 2.3, 1.8, 1.5, 1.7, 3.6, 2.4, 0.0, -1.7, -4.3, -~
$ 기압 <dbl> 1001.3, 1001.9, 1002.6, 1002.8, 1003.0, 1001.8, 1002~
$ 해면기압 <dbl> 1004.9, 1005.5, 1006.2, 1006.4, 1006.6, 1005.4, 1006.3~
$ 일조 <dbl> NA, NA, NA, NA, NA, NA, NA, NA, 0.0, 0.8, 1.0, 1.0, ~
$ 일사 <dbl> 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.05~
$ 전운량 <int> NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, N~
$ 시정 <dbl> 922, 4315, 2601, 1717, 1957, 571, 57, 436, 593, 593,~
$ 지면온도 <dbl> 3.1, 2.8, 2.6, 2.0, 2.2, 4.1, 3.5, 2.2, 1.0, 0.3, 0.6,~
$ floating <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 36, 313, 532, 607, 614,~
$ warehouse <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 33, 209, 296, 315, 474,~
$ dangjin <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 37, 318, 490, 550, 727,~
$ ulsan <int> 0, 0, 0, 0, 0, 0, 0, 0, 4, 35, 71, 82, 334, 372, 3~
time 기온 강수 풍속 풍향 습도 증기압 이슬점 기압
1 0 0 0 0 0 0 0.001835292 0.00191339 0.001718146
해면기압 일조 일사 전운량 시정 지면온도 floating
1 0.001679097 0.4545668 0 0.1552579 0 0 0
warehouse dangjin ulsan
1 0 0 0
울산 지역의 전력 발전량 데이터만 활용할 것이기 때문에 날짜 변수를 제외한 나머지 변수는 제거했다.
ulsan <- rdata %>%
select(-c(dangjin, warehouse, floating)) %>%
select(time, ulsan) %>%
rename(date = time, value = ulsan)
tidymodels의 경우 train/test 분리를 위해서 initial_split()을 활용했는데 시계열 데이터의 경우 특정 날짜를 기준으로 잘라야하기 때문에 Modeltime 패키지에 내장되어있는 initial_time_split()을 이용한다. 특정 날짜를 기준으로 자르고 싶을 경우 timeseries_split()을 이용할 수도 있다.
tk_time_series_cv_plan() : split된 object를 데이터프레임으로 전환
plot_time_series_cv_plan() : sampling된 데이터를 이용해서 시계열 그래프 생성
initial_time_split(data = ulsan, prop = 0.9) %>%
tk_time_series_cv_plan() %>%
plot_time_series_cv_plan(date, value,
.interact = FALSE,
.title = "Partition Train / Test")

months <- 1
total_months <- lubridate::interval(base::min(ulsan$date),
base::max(ulsan$date)) %/%
base::months(1)
prop <- (total_months - months) / total_months
splits <- rsample::initial_time_split(ulsan, prop = prop)
splits %>%
timetk::tk_time_series_cv_plan() %>%
timetk::plot_time_series_cv_plan(date, value)
모델 fitting은 tidymodels 패키지의 방식과 동일하다. 현재는 default 세팅으로 모델을 fitting했지만 튜닝 파라미터가 있을 경우 이전에 tidymodels에서 했던 방식 그대로 grid_latin_hypercube(), grid_random, bayes tuning 등을 이용해서 최적의 파라미터를 찾고 모델을 fitting 해야한다.
# Exponential smoothing
model_fit_ets <- modeltime::exp_smoothing() %>%
parsnip::set_engine(engine = "ets") %>%
parsnip::fit(value ~ date, data = training(splits))
# ARIMA
model_fit_arima <- modeltime::arima_reg() %>%
parsnip::set_engine("auto_arima") %>%
parsnip::fit(
value ~ date,
data = training(splits))
# ARIMA Boost
model_fit_arima_boost <- modeltime::arima_boost() %>%
parsnip::set_engine("auto_arima_xgboost") %>%
parsnip::fit(
value ~ date + as.numeric(date) + month(date),
data = training(splits))
# Prophet
model_fit_prophet <- modeltime::prophet_reg() %>%
parsnip::set_engine("prophet") %>%
parsnip::fit(
value ~ date,
data = training(splits))
# Prophet Boost
model_fit_prophet_boost <- modeltime::prophet_boost() %>%
parsnip::set_engine("prophet_xgboost") %>%
parsnip::fit(
value ~ date + as.numeric(date) + month(date),
data = training(splits))
modeltime_table에 fitting한 모델을 추가한다. modeltime_table은 이전에 서술했다시피 각 모델이 재대로 적합되었는지 확인하고, 이후 예측 워크플로우를 위해서 modeltime_table 구조를 이용하므로 model table을 에러 없이 세팅하는 것이 중요하다.
model_tbl <- modeltime::modeltime_table(
model_fit_ets,
model_fit_arima,
model_fit_arima_boost,
model_fit_prophet,
model_fit_prophet_boost)
model_tbl
# Modeltime Table
# A tibble: 5 x 3
.model_id .model .model_desc
<int> <list> <chr>
1 1 <fit[+]> ETS(A,N,A)
2 2 <fit[+]> ARIMA(5,0,0)(2,1,0)[24]
3 3 <fit[+]> ARIMA(2,0,1)(2,1,0)[24] W/ XGBOOST ERRORS
4 4 <fit[+]> PROPHET
5 5 <fit[+]> PROPHET W/ XGBOOST ERRORS
이전에 만든 modeltime_table을 test 데이터에 적합시켜서 보정을 하는 단계이다.
calibration_tbl <- model_tbl %>%
modeltime::modeltime_calibrate(testing(splits))
calibration_tbl %>%
modeltime::modeltime_accuracy() %>%
flextable::flextable() %>%
flextable::bold(part = "header") %>%
flextable::bg(bg = "#D3D3D3", part = "header") %>%
flextable::autofit()
.model_id | .model_desc | .type | mae | mape | mase | smape | rmse | rsq |
1 | ETS(A,N,A) | Test | 47.69905 | Inf | 2.191688 | 146.0680 | 73.01855 | 0.7655620 |
2 | ARIMA(5,0,0)(2,1,0)[24] | Test | 22.46559 | Inf | 1.032255 | 140.5838 | 46.02106 | 0.7889317 |
3 | ARIMA(2,0,1)(2,1,0)[24] W/ XGBOOST ERRORS | Test | 42.10938 | Inf | 1.934853 | 142.5776 | 58.27198 | 0.7888361 |
4 | PROPHET | Test | 30.82420 | Inf | 1.416318 | 144.1704 | 48.15598 | 0.7779902 |
5 | PROPHET W/ XGBOOST ERRORS | Test | 55.50741 | Inf | 2.550469 | 146.1378 | 67.38202 | 0.7794468 |
calibration_tbl %>%
modeltime::modeltime_forecast(new_data = testing(splits),
actual_data = ulsan,
conf_interval = 0.90) %>%
modeltime::plot_modeltime_forecast(.legend_show = TRUE,
.legend_max_width = 20)